In [4]:
%matplotlib inline
from functools import cache
import numpy as np
import matplotlib.pyplot as plt
import itertools
import math
from decimal import Decimal, getcontext
In [5]:
def make_pq(p1, p2, p3):
    @cache
    def p(i, k):
        if k < 0:
            raise ValueError("p called with k < 0")
        if k == 0:
            return [p1, p2, p3][i - 1]
        return p(i, k - 1) * (p(i, k - 1) + 2 * p(i % 3 + 1, k - 1))
    
    @cache
    def q_odd(i, k):
        return q(i, k - 1) * (q(i, k - 1) + 2 * q(i % 3 + 1, k - 1))
    
    @cache
    def z(i, k):
        return p(i, k) * (p(i, k) + 2 * p((i+1) % 3 + 1, k))
    
    @cache
    def q(i, k):
        if k < 0:
            raise ValueError("q called with k < 0")
        if k == 0:
            return p(i, 0) * (p(i, 0) + 2 * p((i + 1) % 3 + 1, 0))
        return (
            q_odd(i, k) * (z(i, k) + z(i % 3 + 1, 1))
            + q_odd(i % 3 + 1, k) * z(i, k)
        )
    
    return (p, q)
In [6]:
p, q = make_pq(.1, .8, .1)
In [7]:
fig, ax = plt.subplots()
for i in [1, 2, 3]:
    data = []
    for k in range(50):
        data.append(p(i, k))
#         if data[-1] > 1:
#             print(i, k)
#             error
#     print(data[30:40])
    ax.plot(data, label=f"$i={i}$")
ax.legend()
plt.show()
In [8]:
fig, ax = plt.subplots()
for i in [1, 2, 3]:
    data = []
    for k in range(50):
        data.append(q(i, k))
    ax.plot(data, label=f"$i={i}$")
ax.legend()
plt.show()
In [10]:
from datetime import datetime
def plot(p1, p2, p3, k, label=None):
    fig_both, ax_both = plt.subplots()
    fig_p, ax_p = plt.subplots()
    fig_q, ax_q = plt.subplots()
    
    p, q = make_pq(p1, p2, p3)
    colors = ["red", "blue", "green", "pink", "aqua", "lime"]
    color_index = 0
    for fun in [p, q]:
        for i in [1, 2, 3]:
            data = []
            for k_val in range(k):
                data.append(fun(i, k_val))
            (ax_p if fun.__name__ == "p" else ax_q).plot(data, label=f"${fun.__name__}_{i}$", color=colors[color_index])
            ax_both.plot(data, label=f"${fun.__name__}_{i}$", color=colors[color_index])
            color_index += 1
    for fig, ax, name in [[fig_both, ax_both, "both"], [fig_p, ax_p, "p"], [fig_q, ax_q, "q"]]:
        ax.legend()
        ax.set_title(f"$p_1^{{(0)}}={p1}$, $p_2^{{(0)}}={p2}$, $p_3^{{(0)}}={p3}$, $k \\in [0, {k}]$" if label is None else label)
        #fig.savefig(f"./de2/{name}_p1_{p1}-p2_{p2}-p3_{p3}-k_{k}_{datetime.now().strftime('%y-%m-%d-%H-%M')}.png")
        fig.show(warn=False)
In [11]:
plot(.4, .2, .4, 50)
plot(.3, .4, .3, 50)
plot(.2, .6, .2, 50)
plot(.1, .8, .1, 50)
plot(.05, .9, .05, 50)
plot(.025, .95, .025, 50)
In [18]:
getcontext().prec = 1500
delta = 10
plot(Decimal(delta) / Decimal(100), Decimal(100 - 2 * delta) / Decimal(100), Decimal(delta) / Decimal(100), 4096)
In [19]:
getcontext().prec = 3000
delta = 10
plot(Decimal(delta) / Decimal(100), Decimal(100 - 2 * delta) / Decimal(100), Decimal(delta) / Decimal(100), 8192)
In [20]:
getcontext().prec = 1500
epsilon = Decimal(1) / Decimal(2 ** 75)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 4096, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-75}$, $k \in [0, 4096]$")
In [21]:
getcontext().prec = 3000
epsilon = Decimal(1) / Decimal(2 ** 80)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-80}$, $k \in [0, 8192]$")
In [22]:
getcontext().prec = 6000
epsilon = Decimal(1) / Decimal(2 ** 80)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-80}$, $k \in [0, 8192]$")
In [23]:
getcontext().prec = 6000
epsilon = Decimal(1) / Decimal(2 ** 75)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-75}$, $k \in [0, 8192]$")
In [24]:
getcontext().prec = 12000
epsilon = Decimal(1) / Decimal(2 ** 75)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-75}$, $k \in [0, 8192]$")
In [25]:
getcontext().prec = 12000
epsilon = Decimal(1) / Decimal(2 ** 15)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-15}$, $k \in [0, 8192]$")
In [26]:
getcontext().prec = 12000
epsilon = Decimal(1) / Decimal(2 ** 18)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-18}$, $k \in [0, 8192]$")
In [27]:
getcontext().prec = 12000
epsilon = Decimal(1) / Decimal(2 ** 20)
plot(epsilon, Decimal(1) - 2*epsilon, epsilon, 8192, label="$p_1^{(0)} = p_3^{(0)} = \\epsilon = 2^{-20}$, $k \in [0, 8192]$")
In [28]:
getcontext().prec = 1000
for ratio in [(4, 10), (3, 10), (2, 10), (1, 10), (5, 100), (25, 1000)]:
    epsilon = Decimal(ratio[0]) / Decimal(ratio[1])
    middle = Decimal(1) - 2 * epsilon
    plot(epsilon, middle, epsilon, k=50)
In [29]:
getcontext().prec = 1000
for ratio in [(4, 10), (3, 10), (2, 10), (1, 10), (5, 100), (25, 1000)]:
    epsilon = Decimal(ratio[0]) / Decimal(ratio[1])
    middle = Decimal(1) - 2 * epsilon
    plot(epsilon, middle, epsilon, k=100)
In [30]:
getcontext().prec = 1000
for ratio in [(4, 10), (3, 10), (2, 10), (1, 10), (5, 100), (25, 1000)]:
    epsilon = Decimal(ratio[0]) / Decimal(ratio[1])
    middle = Decimal(1) - 2 * epsilon
    plot(epsilon, middle, epsilon, k=300)
In [31]:
getcontext().prec = 10000
for ratio in [(4, 10), (3, 10), (2, 10), (1, 10), (5, 100), (25, 1000)]:
    epsilon = Decimal(ratio[0]) / Decimal(ratio[1])
    middle = Decimal(1) - 2 * epsilon
    plot(epsilon, middle, epsilon, k=300)
In [32]:
getcontext().prec = 1000
p1 = Decimal(4) / Decimal(10)
p3 = Decimal(1) / Decimal(10)
p2 = Decimal(1) - p1 - p3
plot(p1, p2, p3, k=50)
In [33]:
getcontext().prec = 1000
p1 = Decimal(10) / Decimal(100)
p3 = Decimal(1) / Decimal(100)
p2 = Decimal(1) - p1 - p3
plot(p1, p2, p3, k=50)
In [35]:
getcontext().prec = 1000
epsilon = Decimal(2) / Decimal(10)
plot(epsilon, 1 - 2 * epsilon, epsilon, k=50)
In [36]:
getcontext().prec = 1000
p1 = Decimal(10) / Decimal(100)
p3 = Decimal(1) / Decimal(100)
p2 = Decimal(1) - p1 - p3
plot(p1, p2, p3, k=500)
In [5]:
for prec in [1, 2, 4, 8, 16, 32, 64]:
    print(prec * 1_000)
    getcontext().prec = prec * 1_000
    p1 = Decimal(10) / Decimal(100)
    p3 = Decimal(1) / Decimal(100)
    p2 = Decimal(1) - p1 - p3
    plot(p1, p2, p3, k=500)

1000
2000
4000
8000
16000
32000
64000

<ipython-input-3-d627a13c841f>:5: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig_q, ax_q = plt.subplots()
In [8]:
p1, p2, p3 = .1, .8, .1
p, q = make_pq(p1, p2, p3)
for i in range(10):
    p1, p2, p3 = p(1, i), p(2, i), p(3, i)
    q1, q2, q3 = q(1, i), q(2, i), q(3, i)
    print(i)
    print(f"{p1=} {p2=} {p3=}")
    print(f"{q1=} {q2=} {q3=}")
    print()

0
p1=0.1 p2=0.8 p3=0.1
q1=0.030000000000000006 q2=0.8 q3=0.17000000000000004

1
p1=0.17000000000000004 p2=0.8 p3=0.030000000000000006
q1=0.08216799000000005 q2=0.9120000000000004 q3=0.005832010000000004

2
p1=0.30090000000000006 p2=0.6880000000000002 p3=0.011100000000000006
q1=0.23611145499121605 q2=0.7613652402452181 q3=0.002523304763567201

3
p1=0.5045792100000002 p2=0.4886176000000002 p3=0.006803190000000005
q1=0.5650728557122904 q2=0.4318257554034733 q3=0.0031013888842398923

4
p1=0.7476927443644167 p2=0.24539547577004822 p3=0.0069117798655359075
q1=0.9122555293526421 q2=0.08295397757646429 q3=0.004790493070902393

5
p1=0.9260052734414311 p2=0.06361117854545037 p3=0.01038354801312015
q1=0.9888850512370223 q2=0.002014190283815465 q3=0.00910075847918303

6
p1=0.9752943400071625 p2=0.00536740148913685 p3=0.01933825850370388
q1=0.9813464958493339 q2=0.00019028358868335618 q3=0.01846322056203074

7
p1=0.9616686422358091 p2=0.0002364013937257756 p3=0.038094956370471625
q1=0.962005978095993 q2=1.6647141910773146e-05 q3=0.03797737476220497

8
p1=0.9252612570743184 p2=1.80672871787597e-05 p3=0.0747206756385158
q1=0.9203191270868243 q2=2.498332868829839e-06 q3=0.07967837458055035

9
p1=0.8561418277644417 p2=2.7003262367700373e-06 p3=0.14385547190934733
q1=0.8294636397171021 q2=7.157078601871746e-07 q3=0.17053564457557613
In [10]:
for k in [10, 15, 20, 25, 30]:
    plot(.1, .8, .1, k=k)